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Abstract 

We  report  on  the  research  activity  developed  under  the  present  contract 
in  the  field  of  anisotropic  mesh  adaption  and  error  estimation.  With  respect 
to  the  former  problem,  we  describe  a  method  for  the  insertion  of  anisotropic 
mesh  layers  in  general  tetrahedral  grids.  For  a-posteriori  error  estimation,  we 
describe  an  implementation  of  a  recovery-based  estimator  and  its  extention 
to  anisotropic  metric-based  adaption  strategies.  Selected  examples  illustrate 
the  characteristics  of  the  proposed  procedures.  We  conclude  by  commenting 
on  the  results  of  the  research  activity  and  by  discussing  possible  future  work 
in  the  same  field. 


Key  Words 

Finite  elements,  adaptivity,  anisotropic  meshes,  unstructured  tetrahedral 
and  hybrid  grids,  a-posteriori  error  estimation,  recovery-based  error  estima¬ 
tion,  metric-based  grid  adaption. 
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1.  Introduction 

In  this  report  we  describe  the  work  we  have  performed  with  the  support 
of  the  U.S.  Army  Research  Office  through  its  European  Office  in  London 
under  a  contract  with  the  Politecnico  di  Milano.  The  work  here  described 
aims  at  investigating  promising  solution  procedures  for  the  anisotropic  mesh 
adaptive  analysis  of  PDE’s,  with  particular  reference  to  the  compressible 
Euler  and  Navier-Stokes  equations. 

Indeed,  solution  features  in  many  different  branches  of  computational  me¬ 
chanics  are  associated  with  strong  directionality.  In  the  area  of  external  fluid 
dynamics,  for  example,  we  can  cite  the  case  of  rotor  aerodynamics  as  a  typi¬ 
cal  situation  where  the  anisotropy  of  the  solution  plays  a  very  important  role 
and  poses  serious  limitations  to  the  effectiveness  of  isotropic  mesh  methods. 
For  example,  when  refining  shocks  or  pressure  signals  in  high  speed  impul¬ 
sive  noise  analysis,  it  would  be  extremely  desirable  to  have  the  ability  to 
generate  elements  that  are  highly  stretched  in  the  direction  of  the  shock 
surface.  In  the  tracking  of  vortices,  gradients  are  much  higher  in  the  radial 
direction  so  that  highly  stretched  elements  in  the  direction  of  the  vortex 
itself  are  highly  desirable.  On  the  contrary,  isotropic  mesh  refinements  are 
inefficient,  in  the  sense  that  they  tend  to  generate  too  many  grid  points, 
most  of  them  bearing  no  practical  utility.  Since  computer  resources  are 
always  limited  and  since  three  dimensional  problems  are  always  extremely 
demanding  in  terms  of  computing  time  and  memory,  the  lack  of  efficiency 
in  this  sense  seriously  undermines  the  ability  of  such  isotropic  processes  to 
demonstrate  grid  independent  results. 

The  work  is  here  organized  in  two  main  sections.  Section  2  discusses  the 
mesh  modification  procedures  that  are  used  for  inserting  layers  of  anisotropic 
elements  into  three-dimensional  grids.  The  methodologies  here  illustrated 
produce  anisotropic  grids  of  excellent  quality  by  carefully  avoiding  the  con¬ 
struction  of  elements  with  large  dihedral  angles.  This  result  is  achieved  by 
inserting  structured  layers  of  elements  that  are  finally  tetrahedronized. 

Next,  Section  3  discusses  the  problem  of  error  estimation  for  driving  the 
adaption  process.  Two  procedures  are  described.  The  first  makes  use  of 
a  generalization  of  the  Zienkeiwicz-Zhu  (ZZ)  error  estimator,  that  produces 
robust  estimates  at  low  computational  cost.  The  method  is  essentially  based 
on  the  computation  of  higher  order  information  regarding  the  gradients  of 
the  solution.  The  estimate  is  then  obtained  by  comparing  the  original  com¬ 
puted  gradient  and  the  reconstructed  one.  This  procedure  is  known  to  pro¬ 
duce  excellent  results  for  a  number  of  different  applications,  and  it  is  here 
implemented  in  our  research  parallel  adaptive  CFD  code.  One  drawback  of 
this  approach  is  that  no  anisotropic  information  is  directly  available  from  the 
computed  estimator  for  driving  a  metric-based  anisotropic  adaption  process. 
To  circumvent  this  limitation,  a  second  procedure  is  proposed  that  is  based 
on  estimates  of  the  interpolation  error  which  explicitely  take  into  account 
the  anisotropy  of  the  solution  and  of  the  mesh  elements.  These  estimates, 
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which  in  principle  would  be  prohibitively  expensive  to  evaluate  for  realistic 
problems,  are  then  made  computable  by  using  the  Zienkeiwicz-Zhu  recon¬ 
structed  gradients.  This  way  a  local  metric  can  be  constructed  for  driving 
an  anisotropic  iterative  adaptive  mesh  method. 

Finally,  Section  4  presents  some  final*  remarks  and  comments  on  the  work 
conducted  so  far.  In  particular,  we  mention  the  highlights  and  the  present 
limitations  of  the  layer  insertion  algorithms  reported.  We  propose  to  comple¬ 
ment  the  present  procedures  with  metric-based  point  insertion  algorithms, 
that,  at  the  price  of  lower  quality  grids,  will  extend  the  applicability  of  the 
adaptive  software  to  handle  more  general  situations.  Next,  we  briefly  discuss 
the  error  estimation  procedures  and  we  propose  plans  of  future  work  in  this 
area.  Finally,  we  conclude  listing  the  scientific  publications  produced  under 
the  support  of  this  contract. 
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2.  Procedures  for  Anisotropic  Mesh  Adaption 

2.1.  Introduction  and  Motivation.  Many  problems  in  fluid  and  solid 
mechanics  present  solution  features  that  are  inherently  directional.  For  ex¬ 
ample  shock  waves  are  characterized  by  gradients  that  are  very  high  in  one 
local  direction,  while  substantially  smaller  in  the  other  two.  Similarly,  a 
vortex  shed  by  the  tip  of  a  wing  or  a  helicopter  blade  presents  gradients  in 
the  radial  direction  that  are  much  higher  than  in  the  direction  of  its  core. 
In  all  these  and  many  other  cases,  anisotropic  meshes  are  highly  desirable 
for  their  evident  economic  advantages. 

It  appears  that  structured  procedures  are  better  suited  at  generating 
highly  anisotropic  grids.  However,  unstructured  grids  can  discretize  non¬ 
manifold  domains  of  arbitrary  topological  complexity  with  very  little  user 
intervention  and  consequent  limited  work  load.  Furthermore,  unstructured 
grids  can  be  locally  adapted  in  order  to  capture  the  relevant  solution  fea¬ 
tures  to  the  required  level  of  approximation.  These  two  characteristics  of 
the  unstructured  approach  are  contributing  to  its  ever-increasing  success  in 
complex  applications  in  many  diverse  disciplines. 

In  general,  the  introduction  of  anisotropic  regions  can  be  obtained  during 
unstructured  mesh  generation  with  the  use  of  local  mappings  [5],  although 
the  amount  of  element  stretching  is  somewhat  limited  by  the  need  to  pro¬ 
vide  some  form  of  smooth  transition  to  isotropic  regions.  Furthermore,  the 
method  tends  to  create  large  dihedral  angles  that  affect  the  overall  quality 
of  the  resulting  mesh. 

One  possible  alternative  is  represented  by  the  introduction  of  structured 
mesh  regions  in  the  areas  of  localization,  while  using  an  unstructured  dis¬ 
cretization  to  fill  the  rest  of  the  computational  domain.  The  nature  of 
structured  discretizations  makes  it  possible  for  the  contrasting  requirements 
of  very  high  elongation  and  good  mesh  quality,  in  particular  control  of  large 
dihedral  angles,  to  be  more  easily  achieved.  For  example,  ref.  [28]  develops  a 
semi-structured  two-dimensional  procedure  that  first  forms  layers  of  quadri¬ 
lateral  elements  in  the  regions  of  high  localization  and  then  triangulates  the 
rest  of  the  domain  using  an  advancing  front  method.  However,  in  three 
spatial  dimensions  the  use  of  local  structured  hexahedral  meshes  would  lead 
to  non-conformity  at  the  structured/unstructured  interfaces,  with  obvious 
complications. 

In  three  spatial  dimensions,  hybrid  tetrahedral/prismatic  grids  have  been 
used  with  success  in  regions  of  high  gradients  in  proximity  of  model  faces 
(boundary  layers).  The  prismatic  regions  present  a  structured  discretiza¬ 
tion  in  the  direction  of  the  solution  gradients,  which  allows  for  extremely 
high  elongations  with  good  control  of  dihedral  angles.  At  the  same  time, 
the  grid  furnishes  an  unstructured  triangular  discretization  of  the  bound¬ 
ary  layer  face  that  allows  for  a  conforming  interface  with  the  rest  of  the 
isotropic  grid.  Refs.  [17,  23]  develop  a  method  for  the  generation  of  hybrid 
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tetrahedral/prismatic  grids,  where  a  prismatic  mesh  is  obtained  by  march¬ 
ing  away  from  the  model  boundary  an  initial  triangular  surface  mesh.  An 
octree-advancing  front  method  is  then  used  for  filling  the  rest  of  the  domain 
with  tetrahedra. 

The  advancing  layers  method  [21,  10],  in  its  various  flavors,  uses  similar 
ideas.  The  method  starts  with  a  triangulation  of  the  model  faces,  that  is 
then  inflated  along  quasi-normal  directions  using  a  modified  advancing  front 
algorithm.  Various  forms  of  pruning  and  collapsing  of  elements  in  problem 
regions  is  used  for  improving  the  robustness  and  generality  of  the  method. 
The  process  is  terminated  by  filling  the  rest  of  the  domain  with  an  isotropic 
mesh  obtained  via  a  standard  advancing  front  method.  A  generalization  of 
the  method  based  on  multiple  growth  directions  is  detailed  in  ref.  [15].  Since 
multiple  growth  curves  can  be  used  at  any  given  point,  gaps  are  created  in 
the  boundary  layer  mesh.  Special  blends  and  transition  elements  are  imple¬ 
mented  to  hide  the  gaps  generated  by  the  multiple  growth  directions  that 
expose  highly  anisotropic  faces  to  the  mesh  generator.  Both  the  advancing 
layers  and  the  hybrid  tetrahedral/prismatic  method  are  designed  primarily 
for  generating  boundary  layer  meshes;  the  anisotropic  refinement  of  internal 
features  such  as  shock  waves  can  be  performed  by  adding  fictitious  entities 
to  the  underlying  solid  model,  for  example  by  inserting  an  internal  model 
face  to  represent  the  shock  wave  surface. 

In  this  section  of  the  report,  we  present  a  set  of  automated  procedures  for 
the  directional  refinement  of  three-dimensional  unstructured  grids  defined  in 
complex  domains,  with  the  specific  goal  of  allowing  large  values  of  element 
stretching  and  good  control  of  the  mesh  quality.  The  method  is  suitable  both 
for  the  generation  of  boundary  layer  meshes  and  for  the  adaptive  refinement 
of  internal  features  with  local  one-dimensional  anisotropy,  i.e.  with  large 
gradients  in  one  single  local  direction  as  in  the  typical  case  of  shock  waves. 
The  proposed  procedures  start  with  an  initial  isotropic  or  anisotropic  grid 
and  insert  layers  of  prisms  in  the  regions  of  localization.  The  location  of 
the  regions  that  need  to  be  adapted  and  the  stretching  of  the  inserted  layers 
are  assumed  to  be  known  from  the  use  of  an  appropriate  error  estimator  or 
from  user-defined  input. 

For  boundary  layers,  the  procedure  is  normally  driven  by  user-supplied  in¬ 
put.  Attributes  specifying  (possibly  non-uniform)  boundary  layer  thickness, 
number  of  layers  and  their  geometric  distribution  are  associated  with  faces 
in  the  solid  model  wherever  a  boundary  layer  needs  to  be  grown.  The  mesh¬ 
ing  process  starts  by  locally  deflating  an  initial  isotropic  grid,  obtained  with 
any  suitable  mesh  generator.  The  deflation  of  the  mesh  is  obtained  with  a 
specially  designed  mesh  motion  algorithm  based  on  force  control  whose  goal 
is  to  make  room  for  the  boundary  layer  region.  The  highly  stretched  tetra¬ 
hedra  of  the  boundary  layer  region  are  then  created  by  first  inserting  stacks 
of  prisms  in  the  resulting  void,  between  the  initial  and  the  final  positions  of 
the  external  triangulation  of  the  grid,  and  then  by  tetrahedronization  of  each 
prism.  A  simple  algorithm  ensures  the  determination  of  the  correct  prism 
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splitting  template  that  needs  to  be  used  in  each  case  in  order  to  obtain  a 
conforming  grid. 

For  features  within  the  problem  domain  (e.g.  shocks),  the  tetrahedral 
mesh  is  first  cut  along  the  feature.  This  creates  a  triangulated  surface  within 
the  mesh.  Next,  the  two  sides  of  the  cut  are  marched  away  from  each  other 
in  order  to  make  room  for  the  prismatic  mesh.  Once  again,  the  required 
space  is  obtained  by  deformation  of  the  tetrahedral  grid. 

In  both  cases,  the  procedures  are  based  on  the  idea  of  opening  a  gap  in 
the  mesh  that  has  the  same  triangular  discretization  on  the  two  sides,  and 
filling  the  resulting  void  with  stacks  of  prisms.  We  refer  to  this  procedure  as 
“wedging”.  The  meshing  processes  are  based  on  a  solid-model  description 
of  the  domain  which  allows  for  a  consistent  handling  of  the  model  geometry 
and  avoids  precision  losses  in  the  presence  of  curved  boundaries. 

The  proposed  method  achieves  two  main  goals:  first,  it  is  able  to  introduce 
“structured”,  and  therefore  well-behaved,  regions  of  elements  within  arbi¬ 
trary  unstructured  grids;  second,  this  result  is  obtained  with  a  remarkably 
simple  process.  In  fact,  all  advancing  layers  methods  grow  the  boundary 
layer  mesh  in  a  void.  This  is  a  difficult  task  for  highly  anisotropic  meshes, 
and  requires  a  number  of  accompanying  complex  correction  procedures.  On 
the  contrary,  the  proposed  method  makes  careful  use  of  the  deformations 
imposed  on  an  already  exiting  isotropic  grid  together  with  appropriate  mesh 
smoothing  processes  in  order  to  obtain  a  natural  way  of  controlling  the  mesh 
quality  by  automatically  adjusting  the  anisotropic  layer  thickness  and  the 
growth  directions.  Once  calibrated,  the  procedures  are  remarkably  simple, 
being  essentially  based  on  a  robust  mesh  smoother  and  on  a  set  of  prism 
splitting  templates,  in  addition  to  a  few  support  routines.  As  far  as  the 
generation  of  boundary  layers  is  concerned,  the  procedures  here  proposed 
are  not  as  general  as  the  advancing  layers  method  of  ref.  [15],  although  they 
seem  to  offer  a  simpler  alternative  in  many  practical  cases.  Furthermore, 
they  can  also  be  applied  to  the  refinement  of  internal  features,  which  would 
not  be  easy  to  achieve  with  the  advancing  layers  method. 

This  section  is  organized  as  follows.  Some  basic  terminology  used  through¬ 
out  this  work  is  introduced  in  Section  2.2.  Section  2.3  introduces  the  prob¬ 
lem  by  discussing  the  internal  insertion  of  anisotropic  layers.  The  required 
mesh  cutting  is  explained  in  2.3.1.  Section  2.4  gives  a  general  overview  of 
the  boundary  layer  meshing  algorithm,  providing  details  on  the  prism  stack 
creation  and  on  the  wedge  tetrahedronization.  Next,  we  discuss  the  mesh 
motion  algorithm  in  Section  2.5,  together  with  a  line  identification  scheme 
used  for  ordering  during  smoothing.  Details  on  the  practical  implementa¬ 
tion  of  the  gridding  algorithm  are  given  in  Section  2.6,  where  we  analyze 
the  problem  of  visibility  of  the  deflation  directions,  we  comment  on  several 
heuristic  quality  enhancement  procedures  that  we  have  found  useful  and 
effective,  on  the  treatment  of  transition  zones  at  the  layer  boundaries  and 
on  the  reclassification  of  the  mesh  entities.  The  proposed  anisotropic  mesh 
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generation  algorithm  is  tested  on  a  number  of  practical  examples  in  Section 
2.7.  The  quality  of  the  resulting  meshes  is  analyzed  using  geometric  criteria. 

2.2.  Solid  Modeling  and  Terminology.  Any  given  geometric  domain  can 
be  described  by  a  boundary  representation  (BRep)  of  its  topology.  The  ba¬ 
sic  topological  information  is  represented  by  topological  entities  that  bound 
each  other,  in  the  order  regions,  faces,  edges  and  vertices.  An  attribute 
information  associated  with  each  topological  entity  specifies  its  geometry, 
e.g.  a  surface  associated  with  a  face  or  a  curve  associated  with  an  edge. 
Throughout  this  work  we  use  the  capital  letters  R,  F,  E,  V  to  indicate 
regions,  faces,  edges  and  vertices,  respectively,  in  a  solid  model. 

Grids  can  also  be  represented  using  a  boundary  representation,  although 
not  all  four  entities  are  always  necessarily  present  in  any  data  structure. 
In  fact,  the  sets  of  entities  required  for  representing  a  grid  is  somewhat 
dependent  on  the  application  and  the  use  that  the  application  is  making  of 
the  grid  itself  [4].  In  any  case,  all  topological  entities  in  a  mesh  are  generated 
by  discretizing  a  parent  entity  in  the  geometric  model.  For  example,  a  mesh 
face  can  be  obtained  either  by  discretizing  a  model  face  or  by  discretizing 
a  model  region.  The  term  “classification”  identifies  this  unique  relationship 
between  each  entity  in  the  mesh  with  an  entity  in  the  model.  The  symbols 
r,  f,  e,  v  are  used  to  indicate  regions,  faces,  edges  and  vertices,  respectively, 
in  a  mesh. 

The  procedures  described  in  this  work  make  use  of  a  mesh  data  structure 
that  is  interfaced  with  a  computer  aided  design  (CAD)  system,  that  provides 
geometric  and/or  topological  interrogations  on  the  nature  of  the  domain 
through  a  limited  set  of  functional  operations.  Topological  interrogations 
are  used  for  ensuring  the  compatibility  of  the  mesh  with  the  model,  while 
geometric  interrogations  are  here  used  for  the  proper  placement  of  mesh 
vertices  on  the  model  boundary. 

2.3.  Meshing  of  Internal  Features.  We  start  by  describing  the  wedging 
algorithm  for  the  case  of  internal  features,  which  represents  the  most  com¬ 
plex  case.  The  case  of  boundary  features  will  then  just  use  a  subset  of  the 
described  procedures. 

At  this  point  of  the  discussion,  we  will  simply  assume  that  a  directional  er¬ 
ror  estimator  is  available  for  the  problem  at  hand.  A  possible  simple  method 
of  computing  the  regions  of  localization  and  the  required  local  stretching  is 
discussed  in  ref.  [28].  In  any  case,  the  procedure  should  be  able  to  locally 
define  a  minimum  edge  size  hmin  in  the  local  direction  of  anisotropy  n  and 
a  maximum  edge  size  hmax  in  the  orthogonal  direction. 

As  a  preliminary  step,  passes  of  edge-based  refinement  are  applied  to  the 
mesh  in  order  to  meet  the  requirement  on  the  local  value  of  hmax ■  Next,  the 
requirement  on  hmin  needs  to  be  enforced.  To  this  end,  we  create  a  surface 
that  gives  a  geometric  approximation  of  the  feature  to  be  refined.  The 
geometry  is  obtained  through  the  use  of  the  error  estimator  that  indicates 
the  grid  points  where  high  values  of  stretching  are  needed.  The  targeted 
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grid  points  axe  fitted  on  to  a  surface,  that  in  this  way  approximates  the 
mid-surface  of  the  region  of  high  gradients. 

2.3.1.  Mesh  Cutting  and  Local  Optimization.  Next,  the  mesh  is  cut  along  the 
surface.  Each  tetrahedron  intersected  by  the  surface  is  split  into  tetrahedra 
according  to  the  corresponding  subdivision  pattern.  In  turn,  this  provides 
a  triangular  conforming  discretization  of  the  surface.  The  resulting  surface 
mesh  is  stored  in  a  list  for  subsequent  use  in  the  wedging  process. 

In  general,  cutting  the  mesh  along  the  localization  surface  will  locally  gen¬ 
erate  elements  of  poor  quality  that  will  negatively  effect  the  final  adapted 
grid.  Furthermore,  the  quality  of  the  surface  mesh  that  will  provide  the 
triangular  faces  for  the  stack  of  prisms  introduced  during  wedging,  will  be 
poor.  Therefore,  local  optimization  processes  must  be  applied  for  restoring 
the  mesh  quality  in  the  region  of  the  cut.  The  local  optimization  is  inexpen¬ 
sive,  since  it  is  applied  to  a  surface  mesh  and  to  the  layer  of  elements  that 
share  mesh  entities  with  it. 

Different  implementations  of  the  local  procedures  are  possible.  The  fol¬ 
lowing  sequence  of  operations  is  robust  and  very  effective  in  our  experience. 
First,  the  list  of  newly  generated  mesh  edges  is  scanned  and  the  split  loca¬ 
tion  of  the  parent  edge  is  examined.  If  the  split  location  deviates  from  the 
edge  mid  point  more  than  a  preselected  value,  the  edge  is  collapsed.  Edge 
collapsing  in  three  dimensions  is  exemplified  in  figure  1.  The  target  vertex 
for  the  collapsing  is  always  a  vertex  that  belongs  to  the  newly  generated  sur¬ 
face  mesh,  otherwise  the  surface  mesh  would  be  distorted.  When  both  edge 
vertices  belong  to  the  surface  mesh,  the  best  target  vertex  for  collapsing  the 
edge  is  chosen  based  on  the  quality  of  the  resulting  retriangulations  of  the 
two  possible  polyhedra.  Classification  against  the  underlying  solid  model 
entities  allows  to  ensure  the  topological  validity  of  the  collapsing  [25]. 

To  further  improve  on  the  quality  of  the  mesh  in  the  region  of  the  cut, 
we  apply  a  sequence  of  local  optimizations  in  the  following  order:  local  edge 
swaps  on  the  surface  mesh,  face  and  edge  splits.  The  whole  procedure  is 
repeated  until  the  local  quality  of  the  grid  is  restored. 

We  have  now  obtained  a  triangular  grid  that  discretizes  the  surface  ap¬ 
proximating  the  feature  that  will  be  refined.  The  mesh  entities  of  this  grid 
are  now  duplicated  and  separated,  generating  a  single  layer  of  prisms.  The 
required  space  for  inserting  the  prisms  is  obtained  by  slightly  deforming  the 
rest  of  the  mesh  using  a  vertex  repositioning  algorithm  based  on  a  ficti¬ 
tious  elastic  problem,  described  in  the  following.  The  single  layer  is  then 
split  in  order  to  create  stacks  of  prisms,  that  are  subsequently  subdivided 
into  tetrahedra.  The  whole  procedure  is  exemplified  in  figure  2,  that  shows 
the  generation  of  the  surface  discretization,  its  opening  and  the  subsequent 
filling. 

2.4.  Boundary  Layer  Meshing  Algorithm.  For  boundary  layers,  the 
procedure  is  normally  driven  by  user-supplied  input.  Attributes  specifying 
(possibly  non-uniform)  boundary  layer  thickness  t,  number  of  layers  ni  and 
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their  geometric  distribution  are  associated  with  faces  in  the  solid  model 
wherever  a  boundary  layer  needs  to  be  grown.  The  meshing  process  starts 
by  duplicating  the  mesh  entities  classified  on  the  model  faces  that  have 
boundary  layer  meshing  attributes.  This  creates  two  identical  overlapping 
face  discretizations:  one,  termed  the  “model  face  triangulation” ,  will  not  be 
affected  by  the  meshing  algorithm  and  will  provide  the  surface  mesh  on  the 
model  boundary;  the  other,  termed  the  “transition  triangulation”,  will  be 
repositioned  during  the  gridding  process  and  will  connect  the  initial  isotropic 
grid  to  the  newly  formed  boundary  layer  mesh. 

At  this  step  of  the  procedure,  each  couple  of  associated  mesh  faces  defines 
a  triangular  prism  of  null  thickness.  Each  prism  is  then  subdivided  into 
a  stack  of  ni  prisms.  Finally  each  prism  in  the  stack  is  tetrahedronized, 
using  the  algorithm  described  in  Section  2.4.1.  The  boundary  layer  mesh 
is  then  locally  inflated,  i.e.  the  transition  triangulation  is  “pushed”  away 
from  the  model  face  triangulation,  with  the  goal  of  reaching  the  desired 
boundary  layer  thickness.  The  stacks  of  prisms  follow  the  motion  of  the 
transition  triangulation,  providing  a  through-the-thickness  discretization  of 
the  boundary  layer. 

The  initial  isotropic  grid  is  deformed  using  a  mesh  repositioning  scheme 
in  order  to  accommodate  the  motion  of  the  transition  triangulation.  The 
grid  deforming  algorithm  uses  fictitious  springs  that  connect  the  mesh  enti¬ 
ties  as  described  in  Section  2.5.  The  scheme  ensures  that  interpenetration 
of  neighboring  mesh  regions  will  be  avoided,  and  therefore  the  mesh  validity 
will  be  guaranteed  even  for  large  motions.  Furthermore,  the  scheme  is  based 
on  force  control,  i.e.  suitable  forces  are  applied  at  the  transition  triangula¬ 
tion  vertices  so  that  the  requested  amount  of  displacement  (layer  thickness) 
is  obtained  at  low-curvature  points.  The  forces  are  applied  along  marching 
directions  that  are  computed  in  order  to  guarantee  the  “visibility”  of  the 
repositioned  vertex  from  all  mesh  faces  in  the  vertex-manifold,  as  detailed 
in  Section  2.6.1. 

Note  that  only  the  vertices  of  the  transition  triangulation  are  forced  to 
move  along  the  prescribed  marching  directions.  All  other  vertices  on  each 
growth  curve,  i.e.  the  lines  that  connect  each  vertex  on  the  model  face 
triangulation  with  its  sibling  on  the  transition  triangulation,  are  free  to 
be  repositioned  during  the  smoothing  steps  in  the  boundary  layer  mesh. 
This  mesh  relaxation,  once  again  due  to  the  presence  of  fictitious  springs 
that  connect  entities  in  the  grid,  provides  an  automatic  way  of  curving  the 
growth  directions  in  regions  with  closely  spaced  faces  or  of  high  curvature. 
Figure  3  shows  the  automatic  curving  of  growth  directions  obtained  in  the 
common  case  of  a  corner  region. 

The  procedure  automatically  deals  with  problem  configurations,  as  for 
example  concavities  with  high  curvature,  and  avoids  self-intersections  when 
layers  are  grown  on  closely  spaced  model  faces.  In  fact  at  problem  locations, 
as  for  example  points  of  high  curvature  or  sharp  corners,  the  applied  force 


13 


will  automatically  generate  smaller  grid  deformations  than  in  “flat”  loca¬ 
tions,  since  the  fictitious  structure  is  locally  stiffer.  Applied  forces  are  also 
adjusted  by  the  algorithm  according  to  the  local  characteristic  dimensions 
of  mesh  regions,  in  order  to  reduce  the  final  displacement  in  areas  where 
elements  are  smaller,  typically  where  curvature  based  refinement  has  been 
used  by  the  isotropic  mesh  generator.  Finally,  in  order  to  obtain  a  smoother 
variation  of  vertex  displacements,  applied  loads  are  corrected  through  a  dis¬ 
tance  weighted  averaging  among  all  vertices  in  a  vertex-manifold.  A  typical 
grid  resulting  from  this  process  is  exemplified  in  figure  4. 

In  the  case  of  the  advancing  layers  method,  self-intersections  are  usually 
created  in  regions  of  high  curvature  or  at  closely  spaced  faces.  Fixing  of 
boundary  layer  intersections  must  then  be  accomplished  by  shrinking  and 
pruning  of  layers.  Furthermore,  a  localization  structure,  e.g.  an  octree,  is 
needed  in  order  to  identify  all  the  self-intersecting  elements  avoiding  0(n2) 
operations.  We  note  that  the  need  for  intersection  checks  and  fixes  is  com¬ 
pletely  avoided  with  the  present  method,  since  the  possibility  of  creating 
intersections  is  avoided  a-priori. 

The  repositioning  of  the  transition  triangulation  is  usually  performed  in 
incremental  steps,  i.e.  a  first  value  of  the  load  is  imposed  on  the  triangu¬ 
lation  vertices  and  the  isotropic  and  boundary  layer  grids  are  repositioned 
using  a  few  smoothing  steps  of  the  mesh  motion  algorithm;  next,  a  second 
load  increment  is  applied  on  the  vertices  followed  by  smoothing.  The  pro¬ 
cess  continues  until  the  desired  layer  thickness  is  obtained.  Throughout  this 
process,  checks  are  performed  in  order  to  ensure  that  the  computed  dis¬ 
placements  do  not  cause  a  vertex  to  cross  the  planes  identified  by  each  mesh 
face  of  the  polyhedral  cavity  containing  it.  If  a  plane  crossing  is  identified, 
the  applied  load  is  decreased  accordingly.  Consequently,  it  is  possible  that 
a  vertex  does  not  reach  the  desired  final  displacement,  stopped  in  its  move¬ 
ment  by  a  neighboring  face.  In  order  to  avoid  local  abrupt  changes  in  the 
displacements  of  the  transition  triangulation  with  consequent  possible  sever 
deformations  of  the  mesh,  which  of  course  have  a  negative  impact  on  the 
final  result,  it  is  beneficial  to  propagate  these  load  changes  to  the  applied 
loads  of  the  neighboring  vertices.  This  is  here  accomplished  by  a  distance 
weighted  averaging  for  all  vertices  in  the  vertex-manifold. 

2.4.1.  Tetrahedronization.  Prisms  in  the  anisotropic  region  are  tetrahedron- 
ized.  In  fact,  in  general  a  hybrid  mesh  that  consists  only  of  tetrahedra  in 
the  isotropic  part  of  the  grid  and  of  prisms  in  the  anisotropic  layer,  can¬ 
not  be  guaranteed,  unless  suitable  constraints  on  the  definition  of  the  layer 
growth  attributes  are  specified,  as  detailed  in  Section  2.6.3.  Clearly,  the 
tetrahedronization  has  to  be  conforming,  in  the  sense  that  prisms  sharing 
a  quadrilateral  face  need  to  use  the  same  of  the  two  possible  splittings  of 
that  face.  All  possible  tetrahedronization  templates  for  a  prism  are  given 
in  figure  5.  Each  template  can  be  symbolically  represented  by  orienting  the 
edges  of  one  of  the  triangular  faces  of  the  prism,  according  to  a  conventional 
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rule,  as  shown  in  the  figure.  Note  that  no  template  induces  a  closed  loop 
orientation  of  the  edges. 

This  fact  can  be  used  for  constructing  a  simple  algorithm  that  ensures 
a  correct  tetrahedronization.  To  this  end,  the  list  of  vertices  in  the  model 
face  triangulation  is  traversed.  Each  visited  vertex  orients  all  edges  that  use 
it  and  that  have  not  yet  been  oriented.  The  orientation  is  always  defined 
according  to  a  pre-specified  rule,  e.g.  in  the  direction  pointing  towards  the 
vertex  being  examined.  The  procedure  is  exemplified  in  figure  6.  Once  all 
vertices  in  the  surface  mesh  have  been  visited,  all  edges  have  been  oriented. 
It  is  easily  shown  that  no  closed  loop  can  be  generated  by  such  process. 
Therefore,  the  computed  orientations  can  be  used  for  determining  the  correct 
prism  splitting  template. 

2.5.  Mesh  Motion  Algorithm.  The  mesh  motion  algorithm  is  used  for 
deflating  the  isotropic  mesh.  A  suitable  design  of  such  algorithm  is  criti¬ 
cal  for  the  performance  of  the  proposed  method.  The  vertex  repositioning 
problem  must  be  solved  within  the  following  constraints.  First,  the  grid 
must  not  be  severely  deformed  in  the  proximity  of  the  anisotropic  layer  in 
order  not  to  degrade  its  quality.  Second,  it  must  provide  an  automatic  way 
of  dealing  with  problem  regions,  as  mentioned  earlier.  Third,  it  must  be 
robust  and  efficient,  and  in  particular  it  must  not  generate  invalid  elements 
even  for  large  amplitude  motions.  In  order  to  satisfy  the  constraints,  the 
following  procedure  was  devised. 

The  algorithm  is  based  on  a  fictitious  elasticity  problem,  based  on  replac¬ 
ing  the  mesh  edges  with  springs.  Following  ref.  [3],  the  spring  stiffnesses 
are  inversely  proportional  to  the  edge  lengths  in  the  initial  configuration, 
so  that  small  elements  tend  to  be  very  stiff  while  larger  elements  are  softer. 
However,  the  method  in  its  classical  form  can  entangle  the  mesh  when  a  ver¬ 
tex  crosses  a  neighboring  mesh  face.  This  means  that  the  classical  algorithm 
can  only  be  used  for  moderate  mesh  displacements.  An  improved  method 
based  on  the  use  of  additional  torsional  springs  at  the  vertices  was  proposed 
in  ref.  [12]  for  two-dimensional  grids.  The  torsional  springs  are  designed 
so  as  to  guarantee  that  neighboring  triangles  can  not  interpenetrate  each 
other,  in  order  to  ensure  that  mesh  entanglement  will  be  avoided.  However, 
the  introduction  of  torsional  springs  renders  the  problem  non-linear,  so  that 
a  linearized  version  of  the  scheme  must  be  developed.  Furthermore,  the 
method  becomes  quite  cumbersome  in  three  spatial  dimensions. 

Hence,  we  propose  here  a  new  method  that  is  simple,  yet  guarantees  valid 
meshes  even  for  motions  of  large  amplitude.  The  basic  idea  is  to  use  the 
classical  edge  springs  of  ref.  [3]  together  with  additional  linear  springs  that 
prevent  a  mesh  vertex  from  crossing  a  neighboring  face.  For  each  region  r 
using  a  vertex  v,  a  new  linear  spring  is  constructed  that  connects  v  with 
its  projection  v'  on  the  plane  of  the  face  f  of  r  opposite  v.  This  additional 
spring  is  exemplified  in  figure  7  for  a  single  tetrahedron  connected  to  a 
vertex,  for  clarity. 
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Once  the  additional  springs  are  constructed  for  all  tetrahedra  using  a 
vertex,  one  has  effectively  constrained  the  same  vertex  not  to  leave  the 
polyhedral  ball  that  encloses  it,  therefore  ensuring  a  final  valid  mesh.  Fur¬ 
thermore,  the  presence  of  the  additional  springs  is  also  beneficial  in  terms  of 
mesh  quality,  since  they  will  tend  to  keep  each  vertex  closer  to  the  centroid 
of  the  ball,  pushing  it  away  from  its  boundary.  The  resulting  mesh  repo¬ 
sitioning  scheme  is  linear,  and  can  be  straightforwardly  modified  to  handle 
the  two-dimensional  case. 

In  classical  mesh  motion  problems,  the  displacements  at  an  interface,  e.g. 
a  fluid/solid  boundary,  are  imposed  at  the  grid  vertices  at  that  interface. 
The  displacements  are  then  propagated  to  the  rest  of  the  grid  using  some 
form  of  smoothing.  This  ensures  that  the  required  displacement  is  always 
achieved,  at  least  as  long  as  the  smoothing  algorithm  does  not  entangle  the 
mesh.  Such  an  approach  is  unsuitable  for  the  present  application.  A  better 
way  of  deforming  the  mesh  is  through  an  algorithm  based  on  force  control, 
not  displacement  control.  In  practice,  a  force  is  applied  at  mesh  vertices 
on  the  transition  triangulation  along  a  suitable  local  direction.  This  force 
is  calibrated  in  such  a  way  as  to  produce  the  required  local  displacement 
in  a  nearly  flat  boundary  case.  When  the  boundary  presents  sharp  corners 
or  high  curvature,  the  local  forces  are  not  able  to  deform  the  mesh  to  the 
extent  of  producing  the  required  displacement,  since  high  local  deformations 
would  be  required.  These  deformations  are  resisted  by  the  stretching  and 
compressing  of  the  linear  springs. 

Care  is  exercised  in  order  to  ensure  the  visibility  of  the  displaced  vertices 
from  their  parents,  and  for  guaranteeing  that  the  deflation  process  advances 
along  directions  close  to  the  local  normal.  The  computation  of  local  growth 
directions  is  discussed  in  Section  2.6.1. 

The  smoothing  process  is  based  on  a  vertex-by-vertex  relaxed  Gauss- 
Seidel  scheme.  Vertices  in  the  isotropic  region  are  ordered  using  the  algo¬ 
rithm  described  in  Section  2.5.1,  with  the  goal  of  maximizing  for  each  newly 
processed  vertex  the  number  of  already  processed  adjacent  vertices.  This 
helps  in  propagating  effectively  the  displacements  throughout  the  mesh  even 
with  a  small  number  of  iterations. 

The  same  smoothing  is  applied  in  the  isotropic  and  anisotropic  parts  of 
the  grid.  However,  while  in  the  isotropic  grid  the  goal  of  smoothing  is 
that  of  deforming  the  mesh  in  order  to  accommodate  the  mesh  layer,  in 
the  anisotropic  mesh  it  determines  the  distribution  law  of  layer  thicknesses 
normal  to  the  wall.  In  fact  the  transition  triangulation,  pushed  by  the  ap¬ 
plied  loads,  pulls  with  it  the  various  layers  of  prisms  away  from  the  model 
boundary,  layers  that,  at  the  beginning  of  the  process,  are  of  null  thickness 
and  axe  all  collapsed  onto  the  model  face  triangulation.  This  pull  is  real¬ 
ized  by  the  smoothing  process  in  the  anisotropic  grid.  A  small  number  of 
iterations  will  produce  thinner  layers  close  to  the  model  face,  approximat¬ 
ing  a  geometric  distribution,  while  a  larger  number  of  iterations  will  tend 
to  create  an  equally  distributed  spacing  between  layers.  This  is  easily  seen 
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in  a  one-dimensional  case,  where  one  has  a  string  of  points  connected  by 
equal  springs.  The  points  are  initially  all  coincident.  Next,  the  first  point 
is  displaced  by  a  distance  d.  All  other  points  are  then  processed  using  a 
Gauss-Seidel  iteration.  The  second  point  in  the  string  is  repositioned,  and 
its  equilibrium  position  is  at  d/2.  Next,  the  third  point  is  considered,  and  its 
equilibrium  position  is  found  at  d/4,  and  so  on.  When  one  full  Gauss-Seidel 
iteration  has  been  performed,  all  points  have  been  processed  once  and  their 
positions  are  given  by  a  geometric  distribution.  Increasing  the  number  of 
Gauss-Seidel  iterations  will  tend  to  equidistribute  the  point  positions,  yield¬ 
ing  equal  spaces  between  them.  It  is  clear  that  other  distributions  could  be 
obtained  using  progressively  stiffer  or  softer,  rather  than  constant,  springs 
connecting  the  points.  This  simple  one-dimensional  model  approximates  the 
behavior  of  the  grid  vertices  along  a  growth  curve  in  the  anisotropic  layer 
mesh.  Note  that  in  this  region,  different  numbers  of  smoothing  steps  can  be 
used  along  growth  curves  or  in  orthogonal  directions.  In  fact,  while  the  for¬ 
mer  controls  the  layer  spacings,  the  latter  curves  the  growth  directions  and 
might  require  a  different  number  of  smoothing  passes  to  yield  a  satisfactory 
solution. 

Using  the  advancing  layers  method,  one  explicitly  constructs  the  layers 
according  to  a  specific  spacing  function  and  would  therefore  probably  not 
use  a  smoother  inside  the  boundary  layer  region.  This  approach  on  the 
one  hand  gives  total  control  on  the  layer  thicknesses,  but  at  the  same  time 
makes  a  little  harder  to  guarantee  the  element  quality  since  there  are  few  de¬ 
grees  of  freedom  to  play  with,  especially  if  one  does  not  use  curved  growth 
directions.  On  the  contrary,  it  is  clear  from  the  previous  discussion  that 
the  proposed  method  can  only  approximate  the  spacing  function  using  the 
smoother.  The  control  on  the  layer  thicknesses  is  then  somewhat  reduced 
in  the  present  case  when  compared  with  the  advancing  layers  technique. 
However,  it  should  be  remarked  at  this  point  of  the  discussion  that  one 
characteristic  of  the  proposed  approach  is  that  all  the  basic  steps  of  the  gen¬ 
eration  of  an  anisotropic  layer  region,  namely  computation  of  the  growth 
directions  and  of  the  growth  displacements,  curving  of  growth  curves,  com¬ 
putation  of  the  grid  spacing  normal  to  the  wall,  etc.,  are  now  all  blended 
in  one  single  process,  the  mesh  smoothing.  This  means  that  the  various 
steps  are  not  realized  in  any  particular  order,  as  with  other  algorithms,  but 
are  all  coupled  together  and  in  this  sense  each  one  of  them  affects  all  the 
others.  Therefore,  while  one  is  pushing  on  the  transition  triangulation  for 
making  room  for  the  anisotropic  layer,  one  is  also  relaxing  the  positions 
of  the  new  grid  points  within  the  layer  according  to  the  fictitious  elasticity 
problem  solved  by  the  smoothing  algorithm,  achieving  at  once  all  the  various 
necessary  goals.  While  on  the  one  hand  this  necessarily  reduces  somewhat 
the  available  control  as  previously  noted,  the  intimate  coupling  of  the  ba¬ 
sic  logical  steps  of  the  mesh  generation  process  implied  by  this  approach 
would  seem  to  offer  the  potential  for  a  beneficial  effect  on  the  quality  of  the 
resulting  grids. 


2.5.1.  Mesh  Ordering  Algorithm.  For  efficiently  propagating  the  information 
during  the  smoothing  process,  vertices  need  to  be  ordered.  The  problem  is 
somewhat  similar  to  the  generation  of  the  line  orderings  used  by  relaxation 
algorithms  for  the  iterative  solution  of  discrete  operators  on  unstructured 
grids.  In  view  of  the  fact  that  a  very  limited  number  of  smoothing  steps  will 
be  performed  for  efficiency  reasons,  we  would  desire  to  maximize  for  each 
vertex  the  number  of  already  processed  neighbors. 

A  simple  ordering  scheme  was  devised  for  this  application,  based  on  a 
modified  greedy  algorithm  [13].  First,  model  faces  with  growth  attributes 
are  collected  in  connected  sets,  i.e.  if  two  model  faces  share  a  common 
model  edge  or  vertex  they  are  assigned  to  the  same  set.  Next,  a  mesh  vertex 
classified  on  model  boundary  is  selected  for  each  connected  set.  This  vertex 
represents  the  “seed”  point  for  an  ordering  process.  The  seed  vertex  is 
assigned  the  number  1  and  labeled  vi.  Starting  from  vi,  all  edge-connected 
vertices  that  are  also  classified  on  the  connected  set  or  its  boundaries  are 
gathered,  numbered  and  stored  in  a  linked  list.  The  vertex  that  has  been 
assigned  the  number  2,  v2,  is  now  considered.  The  process  is  repeated,  and 
all  edge-connected  vertices  to  v2  that  are  also  classified  on  the  connected  set 
or  its  boundaries  are  collected  in  the  list  and  numbered.  At  the  end  of  this 
first  phase,  all  the  numbered  vertices  are  on  the  transition  triangulation,  and 
therefore  are  all  directly  affected  by  the  imposition  of  the  loads  that  drive 
the  smoothing  procedure.  The  process  is  now  restarted  from  vertex  vi. 
All  edge-connected  vertices  that  have  not  yet  been  numbered  are  gathered, 
numbered  and  stored  in  the  linked  list.  Next,  vertex  v2  is  considered,  and  all 
its  edge-connected  vertices  that  have  not  yet  been  numbered  are  processed. 
The  procedure  continues  until  all  vertices  have  been  numbered. 

The  algorithm  essentially  constructs  successive  layers  of  vertices,  the  first 
layer  being  the  transition  triangulation,  that  wrap  the  connected  set  and 
march  away  from  it.  The  procedure  is  illustrated  in  figure  8.  A  typical 
example  of  the  ordering  constructed  in  this  fashion  is  given  in  figure  9  for  the 
two-dimensional  mesh  of  a  submarine  model.  The  growth  set  is  composed  of 
all  the  model  edges  that  form  the  submarine  silhouette,  while  the  seed  point 
is  labeled  P.  Colors  represent  the  numbers  associated  with  the  vertices  by 
the  numbering  process.  Note  how  the  ordering  creates  layers  that  enclose 
the  model  boundary  propagating  outwards  towards  the  rest  of  the  domain. 
The  described  process  is  straightforwardly  generalized  to  handle  the  internal 
layer  case. 


2.5.2.  Improved  Centroidal  Smoothing.  An  additional  way  of  locally  improv¬ 
ing  the  grid  quality  is  provided  by  some  variant  of  the  well  known  Laplacian 
smoothing  process  [14].  A  general  form  of  the  scheme  can  be  written  as 
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where  P  is  the  point  associated  with  a  given  vertex  v  and  the  Pi  s  are  the 
points  associated  with  the  vertices  that  are  edge-connected  to  v.  This  form 
of  local  smoothing  was  used  in  the  present  work  as  a  further  help  in  trying 
to  remove  problematic  configurations. 

This  method  is  quite  effective  for  improving  the  grid  locally.  However, 
some  caxe  should  be  exercised  since  the  polyhedral  region  defined  by  the 
Pi's  might  be  non-convex.  This  situation  is  most  likely  to  happen  in  highly 
distorted  regions,  or  in  the  proximity  of  closely  spaced  faces  and  corners, 
exactly  where  the  need  for  local  improvement  is  the  highest.  Since  even  a 
non  optimal  repositioning  might  lead  to  some  local  improvement  in  these 
pathological  situations,  a  modified  strategy  was  devised.  The  basic  idea  is 
to  “regularize”  the  problem  by  constructing  a  convex  region  around  vertex 
v  using  some  new  points  Pi  instead  of  the  vertex  points  Pi-  At  this  point, 
the  update  defined  by  (1)  can  be  used  safely. 

To  present  the  main  idea  behind  the  proposed  algorithm,  we  can  start 
from  the  simpler  two-dimensional  case.  The  edge-connected  points  Pi  define 
in  this  case  a  polygon.  Each  polygon  edge  divides  the  plane  into  two  semi¬ 
planes,  one  that  contains  v  and  that  therefore  defines  valid  repositionings  of 
the  same  vertex,  and  a  second  one  that  will  produce  invalid  repositionings. 
The  “right”  final  convex  polygon  for  the  safe  application  of  one  pass  of 
Laplacian  smoothing  is  now  given  by  the  boundary  of  the  figure  obtained 
by  intersecting  all  the  valid  semi-planes.  This  is  exemplified  in  figure  10  for 
a  simple  concave  configuration.  In  reality  we  do  not  need  to  define  the  real 
edges  of  this  new  regularized  polygon,  but  only  its  vertices,  since  only  these 
are  needed  for  doing  Laplacian  smoothing. 

The  position  of  these  vertices  can  be  practically  computed  by  considering 
subsequent  intersections  between  couples  of  edges  of  the  original  polygon. 
For  a  convex  polygon,  all  these  intersections  coincide  with  the  vertices  of 
the  polygon  itself.  However,  for  a  non-convex  figure,  additional  points  will 
be  generated.  All  points  need  now  to  be  checked  against  all  polygon  edges 
to  determine  whether  they  are  on  the  valid  semi-plane  or  not.  This  is  easily 
determined  by  defining  an  arbitrary  normal  to  the  edge  and  considering  two 
distance  vectors,  one  between  one  of  the  edge  vertices  and  the  vertex  that 
will  be  repositioned,  v,  and  the  second  between  the  same  edge  vertex  and 
the  point  that  needs  to  be  checked.  One  then  computes  the  dot  products 
of  the  edge  normal  with  each  distance  vector;  if  the  two  products  have  the 
same  sign,  the  checked  point  and  vertex  v  are  on  the  same  semi-plane  and 
the  point  will  be  retained,  otherwise  they  are  on  different  sides  and  the 
point  needs  to  be  rejected.  The  process  is  continued  untill  all  points  have 
been  checked.  At  the  end  of  the  procedure,  the  remaining  points  define  the 
regularized  polygon. 

In  the  three-dimensional  case,  the  same  idea  can  be  generalized  by  con¬ 
sidering  now  all  points  generated  by  the  intersections  of  all  the  edges  of  the 
surface  of  the  polyhedral  region  with  all  its  triangular  faces.  Once  again, 
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each  newly  generated  point  needs  to  be  checked  to  see  whether  it  lies  on  the 
valid  semi-space  with  respect  to  vertex  v. 

2.6.  Implementational  Issues. 

2.6.1.  Computing  Deflation  Directions.  Marching  directions  are  based  on 
the  “normal”  to  the  vertex-manifold,  i.e.  to  the  group  of  mesh  faces  using 
the  vertex  and  classified  on  model  faces  or  on  the  cutting  surface.  In  fact,  in 
order  not  to  create  invalid  elements,  the  new  vertex  position  must  be  visible 
from  all  faces  in  the  manifold.  We  note  that  we  could  use  the  normals  to 
the  model  faces  rather  than  the  mesh  faces  when  growing  boundary  layers, 
since  the  procedures  have  access  to  the  true  problem  geometry  through  a 
solid  modeling  system.  However,  for  coarse  meshes  this  might  violate  the 
visibility  condition.  In  this  work,  we  use  the  method  suggested  in  ref.  [18] 
for  computing  the  initial  values  of  the  directions.  The  algorithm  first  renders 
the  problem  two-dimensional  by  finding  the  bisection  plane  for  the  wedge 
identified  by  the  two  faces  fx  and  f2  in  the  v  vertex-manifold  that  form  the 
smallest  dihedral  angle.  Then,  the  visibility  angle  on  that  plane  is  bisected 
to  identify  the  manifold  normal.  A  visibility  cone  can  now  be  constructed 
at  v,  tangent  to  fi  and  f2. 

In  regions  with  high  curvature,  at  corner  points  or  in  proximity  of  model 
faces  forming  small  dihedral  angles,  it  is  beneficial  to  let  the  marching  di¬ 
rections  slightly  deviate  from  the  manifold  normals  computed  as  explained 
above.  This  reduces  the  convergence  or  divergence  of  the  growth  curves, 
therefore  minimizing  the  local  compressing  and  stretching  of  the  transition 
triangulation.  To  this  end,  the  nominal  normals  are  smoothed  by  distance- 
weighted  averaging.  Given  a  vertex  and  its  normal,  first  all  normals  at  the 
vertices  in  the  vertex-manifold  are  smoothed,  i.e.  a  first  ring  of  normals 
around  the  vertex  is  processed.  Next,  all  faces  adjacent  to  the  vertex- 
manifold  are  considered,  and  their  normals  at  vertices  are  smoothed,  i.e. 
a  second  ring  of  normals  around  the  vertex  is  processed.  The  process  is 
continued  until  a  preselected  number  of  rings  (typically  four  or  five)  has 
been  affected  by  the  smoothing.  After  smoothing  at  each  vertex,  the  new 
marching  direction  is  verified  to  lie  within  the  visibility  cone,  in  order  to 
ensure  the  validity  of  the  subsequent  repositioning. 

We  note  once  again  that  the  deflation  directions  determine  the  marching 
vectors  of  the  sole  vertices  in  the  transition  triangulation.  All  other  vertices 
within  the  anisotropic  layer,  i.e.  all  vertices  along  the  growth  curves,  will  be 
allowed  to  deviate  from  these  directions  according  to  the  layer  smoothing 
process.  This  will  provide  an  automatic  way  of  curving  the  growth  curves 
according  to  the  local  problem  geometry. 

2.6.2.  Mesh  Quality  Enhancements.  Mesh  smoothing  is  in  some  sense  lim¬ 
ited  by  the  topology  of  the  grid,  since  it  only  deals  with  the  repositioning 
of  the  vertices,  and  can  not  change  the  local  connectivity.  In  order  to  relax 
this  constraint  and  to  further  ensure  that  the  isotropic  grid  quality  is  not 
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affected  by  the  repositioning  process,  local  retriangulation  operations  are 
recursively  applied  where  needed.  The  local  optimization  criterion  is  based 
on  the  largest  dihedral  angles  of  each  tetrahedron.  The  retriangulation  does 
not  affect  the  topology  of  the  transition  triangulation,  not  to  interfere  with 
the  highly  stretched  elements  in  the  boundary  layer,  unless  severe  local  de¬ 
formations  can  not  be  corrected  through  the  already  described  repositioning 
procedures. 

Five  basic  retriangulations  are  considered:  edge  removal  and  its  inverse 
operation  multi-face  removal,  edge  swapping,  edge  collapsing,  and  entity 
(region,  face  and  edge)  splittings  [11].  As  a  first  step,  the  dihedral  angles 
of  all  elements  considered  for  optimization  are  calculated.  Each  element 
violating  a  user-defined  threshold  value  is  put  into  a  linear  list.  Depending 
on  the  configuration  of  each  element  in  the  list  (number  of  dihedral  angles 
above  the  threshold  value,  topological  or  geometrical  constraints,  etc.),  a 
suitable  subset  of  the  above  mentioned  optimization  procedures  is  applied 
to  eliminate  that  element  in  favor  of  improved  elements.  The  procedures 
might  fail  for  a  specific  element  if  the  resulting  configuration  is  topologically 
or  geometrically  not  valid,  or  if  they  lead  to  a  degradation  of  the  quality 
of  any  element  involved  in  that  local  retriangulation.  In  this  case,  or  if  the 
element  is  improved  but  the  largest  dihedral  angle  is  still  above  the  threshold 
value,  the  element  is  considered  for  improvement  in  a  second  pass,  after  all 
elements  have  been  processed.  Since  the  neighborhood  of  the  elements  that 
failed  in  the  first  pass  may  have  been  modified,  it  is  possible  that  they  are 
fixed  in  a  second  pass.  The  procedure  is  repeated  until  a  given  threshold 
value  is  reached  or  no  further  improvement  can  be  achieved. 

Local  coarsening  can  also  be  used  in  a  different  context.  In  fact,  if  bound¬ 
ary  layers  need  to  be  grown  on  closely  spaced  model  faces  facing  each  other, 
the  squeezing  and  compressing  of  the  isotropic  grid  could  not  allow  to  reach 
the  desired  layer  thicknesses.  In  this  case,  local  coarsening  can  be  applied 
in  the  isotropic  region,  for  easing  the  growth  process. 

2.6.3.  Ending  Layers.  In  the  case  of  hybrid/prismatic  mesh  generators,  lay¬ 
ers  are  grown  around  all  boundary  faces  of  a  body.  In  many  applications 
however,  boundary  layers  are  needed  only  on  a  limited  set  of  the  model 
faces.  The  procedure  here  described  can  deal  with  situations  involving  adja¬ 
cent  faces  with  different  boundary  layer  growth  attributes.  Figure  11  shows 
the  common  case  of  a  boundary  layer  whose  thickness  goes  to  zero  at  a 
model  edge.  This  is  simply  handled  by  specifying  a  null  growth  for  the 
mesh  vertices  classified  on  the  model  edge,  followed  by  the  already  described 
prism  splitting  and  tetrahedronization.  A  subsequent  pass  of  mesh  collaps¬ 
ing  eliminates  all  edges  of  zero  length,  providing  the  required  grid.  Figure 
12  shows  the  other  common  case  of  a  boundary  layer  ending  at  neighboring 
model  faces.  Interaction  with  the  solid  modeler  ensures  the  correct  posi¬ 
tioning  of  the  boundary  layer  vertices  on  curved  model  faces,  as  detailed  in 
Section  2.6.5.  The  selection  between  an  edge-ending  layer  (figure  11)  or  a 
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face-ending  layer  (figure  12)  is  determined  based  on  the  maximum  allowable 
dihedral  angle  in  the  mesh,  typically  160  deg. 

Clearly,  if  all  model  faces  are  tagged  with  a  boundary  layer,  or  if  all  grown 
layers  end  on  model  faces  as  in  the  case  of  figure  12,  the  final  prism  tetra- 
hedronization  can  be  skipped,  and  one  obtains  a  mixed  hybrid/prismatic 
grid. 

There  are  however  situations  that  require  special  attention.  Figure  14(a) 
gives  an  example  of  a  possible  problematic  configuration:  in  this  case  a 
boundary  layer  will  be  grown  on  face  Pi,  but  not  on  faces  F2,  F3  and  F4. 
One  possible  solution  would  be  to  allow  for  two  growth  curves  at  vertex  Vi, 
one  classified  on  edge  Ei  and  one  on  edge  E2-  A  simpler  solution  is  sketched 
in  figure  14(b)  and  (c).  In  essence,  the  boundary  layer  is  allowed  to  spillover 
to  the  neighboring  faces  in  a  transition  region  of  thickness  comparable  to 
that  of  the  boundary  layer  itself.  Figure  14(b)  shows  the  patch  of  mesh 
faces  involved  in  the  process.  Figure  14(c)  shows  the  detachment  of  the 
patch  of  faces  from  the  model  boundary,  and  the  creation  of  the  boundary 
layer  region.  Now  all  boundary  layer  elements  belong  to  one  of  the  two  cases 
of  figure  11  and  12,  so  that  the  standard  procedures  can  be  used. 

Another  example  of  anisotropic  layer  boundaries  that  need  special  han¬ 
dling  is  represented  by  the  case  of  an  internal  layer  that  ends  within  the 
domain,  i.e.  it  does  not  extend  all  the  way  to  a  model  boundary.  A  common 
such  case  is  encountered,  for  example,  when  refining  a  shock  wave  at  the  tip 
of  a  blade;  the  anisotropic  layer  in  this  case  will  end  at  a  certain  distance 
from  the  blade  surface,  close  to  the  shock  boundaries,  and  will  clearly  not 
have  to  extend  to  the  far  field  for  efficiency.  In  this  case,  the  layer  thickness 
will  have  to  be  null  at  its  internal  boundary,  similarly  to  the  previously  dis¬ 
cussed  boundary  layer  case  with  zero  thickness  at  a  model  edge.  Here  again 
a  null  growth  length  is  specified  for  the  mesh  vertices  at  the  internal  edge 
of  the  layer,  followed  by  the  already  described  prism  splitting  and  tetrahe- 
dronization.  A  subsequent  pass  of  mesh  collapsing  eliminates  all  edges  of 
zero  length,  providing  the  required  grid.  The  result  of  such  procedure  is 
exemplified  in  figure  13. 

2.6.4.  Reclassification.  Certain  applications  will  require  the  entities  of  the 
final  adapted  grid  to  be  completely  classified,  for  example  when  the  anal¬ 
ysis  attributes  (boundary  conditions,  loads,  material  properties,  etc.)  are 
associated  with  entities  in  the  solid  model,  so  that  they  can  be  consistently 
transferred  to  the  new  adapted  grid  for  use  by  the  solver.  Furthermore, 
a  correct  classification  information  is  necessary  for  repositioning  the  mesh 
vertices  on  curved  model  boundaries  in  order  to  provide  geometrically  consis¬ 
tent  grids  and  avoid  loss  of  precision,  as  detailed  in  Section  2.6.5.  Because 
of  these  reasons,  it  is  necessary  to  compute  the  classification  of  all  newly 
created  mesh  entities. 

Vertices  classified  on  a  model  entity  of  degree  k  are  always  reclassified  on 
an  entity  of  degree  k  +  1.  Therefore,  vertices  classified  on  model  faces  are 


22 


always  classified  internal.  For  vertices  classified  on  model  edges,  one  has  to 
check  whether  both  model  faces  that  use  that  edge  are  associated  with  a 
boundary  layer  attribute.  In  case  they  both  have  attributes,  the  vertex  is 
classified  internal,  otherwise  it  is  classified  on  the  model  face  that  does  not 
have  the  boundary  layer  attribute.  This  situation  is  exemplified  in  figure 
15.  Mesh  vi  is  classified  on  model  edge  Ei,  which  is  used  by  face  Fi  and 
face  F3.  Since  a  boundary  layer  is  grown  only  on  F3  and  assuming  a  face¬ 
ending  layer,  the  newly  generated  vertex  vi  is  necessarily  classified  on  Fi. 
Similarly,  for  vertices  classified  on  model  vertices  one  has  to  find  the  model 
edge  that  is  used  by  model  faces  that  are  not  tagged  for  boundary  layer 
growing.  Again  with  reference  to  figure  15,  the  newly  generated  vertex  V2 
is  necessarily  classified  on  model  edge  E2. 

2.6.5.  Dealing  with  Curved  Boundaries.  For  accuracy,  it  is  critical  that  newly 
generated  or  repositioned  mesh  vertices  are  placed  on  the  true  model  bound¬ 
ary.  The  classification  information  together  with  the  direct  interfacing  of  the 
procedures  with  a  solid  modeler,  provide  for  a  straightforward  solution  to 
this  problem.  In  fact,  all  mesh  entities,  including  the  ones  generated  in  the 
boundary  layer,  are  classified  against  the  various  model  entities,  as  discussed 
in  Section  2.6.4.  Then,  after  each  smoothing  step,  each  mesh  vertex  clas¬ 
sified  on  model  face  or  model  edge  is  snapped  to  the  boundary.  The  new 
location  is  provided  by  the  solid  modeler  through  a  closest  point  interroga¬ 
tion.  The  position  of  the  mesh  vertex  before  the  smoothing  pass  is  used  as 
the  initial  guess  for  the  interrogation. 

It  should  however  be  noted  that  the  simple  snapping  of  vertices  to  the 
boundary  can  in  some  cases  generate  invalid  or  severely  distorted  elements. 
More  sophisticated  techniques  can  be  used  for  dealing  with  this  issue,  for 
example  by  applying  local  retriangulations  as  suggested  in  ref.  [20]. 

2.7.  Numerical  Examples.  In  this  section,  two  test  cases  demonstrate  the 
developed  procedures  for  the  refinement  of  boundary  features,  while  a  third 
case  shows  the  insertion  of  internal  anisotropic  layers.  More  examples  and 
additional  details  can  be  found  in  refs.  [6,  7]. 

2.7.1.  Mesh  Quality  Measures.  For  the  examples  here  considered,  mesh  qual¬ 
ity  measures  axe  given  in  terms  of  geometric  criteria.  This  allows  to  deter¬ 
mine  the  grid  characteristics  in  a  solver  and  application-independent  man¬ 
ner,  while  actual  numerical  simulations  would  raise  questions  concerning  the 
particular  finite  element  or  finite  volume  formulation,  boundary  conditions 
used  and  other  problem  dependent  issues. 

Two  quality  measures  are  considered  in  the  following.  The  first  is  the  max¬ 
imum  dihedral  angle  between  two  neighboring  mesh  faces,  Qs  =  max^i^  Si, 
where  the  angle  at  edge  i  shared  by  faces  j  and  A;  is  Si  —  7r±  arccos  (rij  ■  n^) , 
with  ni  the  normal  to  face  l.  Large  dihedral  angles  can  negatively  impact 
the  solution  of  partial  differential  equations,  affecting  the  discretization  error 
and  the  conditioning  of  the  discrete  problem  [2].  The  second  quality  measure 
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is  the  ratio  of  the  radii  of  the  inscribed  and  circumscribed  spheres,  Qs  =  r/ R, 
and  gives  a  measure  of  the  stretching  of  an  element.  The  boundary  layer 
meshing  process  should  not  significantly  affect  both  measures  for  the  in¬ 
coming  isotropic  mesh,  while  it  should  generate  elements  in  the  anisotropic 
regions  that  have  the  largest  dihedral  angle  close  to  ir/2,  and  that  match 
the  requested  value  of  stretching. 

2.7.2.  Nacelle.  The  first  example  deals  with  the  growth  of  a  boundary  layer 
on  an  aeronautical  engine  nacelle.  For  this  problem  10  anisotropic  layers 
were  grown  on  the  nacelle  model  faces,  with  geometric  thickness  distribu¬ 
tion  in  the  direction  normal  to  the  wall.  A  constant  boundary  layer  thickness 
was  requested,  implying  a  value  of  stretching  ranging  from  about  5  to  about 
40.  Stretching  varies  because  of  non-uniform  layer  thickness  and  because 
the  model  face  triangulation  is  graded  due  to  curvature  based  refinement. 
The  initial  isotropic  grid  had  119240  tetrahedra,  and  5318  mesh  faces  clas¬ 
sified  on  model  faces  with  growth  attributes.  The  final  resulting  grid  has 
429020  tetrahedra,  and  is  depicted  in  figure  16.  Elements  in  certain  parts 
of  the  boundary  layer  were  removed  to  expose  the  anisotropic  grid  and  the 
underlying  model  face  triangulation.  Elements  in  the  isotropic  part  of  the 
mesh  are  not  shown  for  clarity. 

Figure  17  shows  a  histogram  of  the  maximum  dihedral  angles  for  the 
initial  and  adapted  grids.  For  the  adapted  grid,  only  dihedral  angles  for 
tetrahedra  in  the  isotropic  region  are  shown.  Therefore,  the  figure  compares 
the  quality  of  the  isotropic  grid  as  produced  by  the  mesh  generator  and 
after  it  has  been  affected  by  the  mesh  repositioning  required  for  inserting 
the  boundary  layer.  It  appears  from  the  plot  that  the  distribution  of  angles 
is  even  better  for  the  adapted  than  for  the  initial  grid. 

Figure  18  shows  a  histogram  of  the  maximum  dihedral  angles  for  the  sole 
anisotropic  part  of  the  adapted  grid.  Note  how  the  maximum  angles  are 
clustered  around  7r/2  in  this  region,  deviations  from  this  value  being  due 
to  the  regions  of  high  curvature  close  to  the  lip  of  the  nacelle.  Figure  19 
shows  a  histogram  of  the  ratio  of  the  radii  of  the  inscribed  and  circumscribed 
spheres,  Qs,  for  the  sole  anisotropic  part  of  the  adapted  grid.  The  stretching 
measures  are  clustered  around  the  range  of  expected  values,  from  5  to  40. 

2.7.3.  Submarine.  The  second  example  shows  the  growth  of  a  boundary 
layer  around  a  more  topologically  complex  submarine  model.  In  this  case  6 
anisotropic  layers  were  grown  on  the  submarine  model  faces,  with  geomet¬ 
ric  thickness  distribution  in  the  direction  normal  to  the  wall.  The  initial 
isotropic  grid  had  165305  tetrahedra,  with  6994  mesh  faces  classified  on 
model  faces  with  growth  attributes.  The  final  grid  has  291197  tetrahedra, 
and  is  shown  in  figure  20.  Even  in  this  case,  certain  parts  of  the  boundary 
layer  were  removed  to  expose  the  anisotropic  grid  and  the  underlying  model 
face  triangulation. 

Figure  21  shows  details  of  the  anisotropic  grid.  Both  the  junction  of 
the  leading  edge  of  the  fin  with  the  submarine  tower  and  the  tip  of  the  fin 
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axe  shown.  Various  features  of  the  proposed  procedures  are  simultaneously 
visible,  in  particular  the  tilting  of  growth  directions,  the  curving  of  growth 
curves  and  the  variable  boundary  layer  thickness. 

Even  in  this  case,  the  quality  of  the  initial  grid  is  not  significantly  affected 
by  the  repositioning  required  for  the  insertion  of  the  boundary  layer.  Figure 
22  shows  a  histogram  of  the  maximum  dihedral  angles  for  the  initial  and 
adapted  grids,  in  this  latter  case  considering  tetrahedra  in  the  isotropic 
region  only. 

2.7.4.  Insertion  of  a  Thin  Layer.  The  third  example  concerns  the  introduc¬ 
tion  of  an  internal  anisotropic  layer.  We  consider  a  simple  cubic  domain 
with  a  quarter  sphere  cut-out  at  one  of  the  corners.  First,  a  boundary  layer 
mesh  is  grown  on  the  sphere  surface,  as  previously  explained.  Next,  we  in¬ 
sert  a  thin  layer  of  elements  in  front  of  the  sphere  using  the  internal  layer 
algorithm  detailed  in  the  previous  pages.  This  example  is  entirely  of  an 
academic  nature,  especially  since  no  particular  attention  was  paid  to  pro¬ 
duce  the  necessary  mesh  gradation  for  a  realistic  computation,  but  might  be 
representative  of  the  refinement  of  a  strong  shock  in  front  of  a  blunt  body. 

Figure  23  shows  a  view  of  the  internal  surface  of  the  sphere,  of  the 
anisotropic  layer  and  of  the  isotropic  mesh.  Elements  in  the  isotropic  and 
anisotropic  grids  were  removed  in  order  to  show  the  mesh  interior.  Figure  24 
offers  a  slightly  different  view  of  the  same  problem,  where  now  only  elements 
in  the  isotropic  grid  were  removed. 

Although  the  problem  topology  is  quite  simple  in  this  case,  this  problem 
shows  some  of  the  typical  difficulties  encountered  with  more  realistic  appli¬ 
cations,  namely  the  arbitrary  cut  of  elements  and  the  collision  of  the  cut 
front  with  the  model  boundaries.  Even  in  this  case,  it  was  observed  that 
the  mesh  quality  measures  remained  essentially  unchanged  by  the  proposed 
procedure  in  the  isotropic  mesh,  while  the  anisotropic  grid  is  of  very  high 
quality  with  virtually  no  large  dihedral  angles. 
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3.  Procedures  for  Anisotropic  Error  Estimation 

3.1.  Introduction  and  Motivation.  The  a-posteriori  estimation  of  the 
error  associated  with  finite  element  solutions  is  still  an  area  of  very  ac¬ 
tive  research.  Among  the  many  procedures  proposed,  we  can  mention  here 
residual-based  methods  and  various  types  of  recovery-based  estimators  (see 
refs.  [27,  1]  for  a  review  of  available  techniques). 

A  procedure  that  is  simple  and  yet  very  effective  is  the  Zienkeiwicz-Zhu 
(ZZ)  error  estimator  [29,  30],  that  reconstructs  smoother  gradients  than  the 
ones  associated  with  the  numerical  solution  of  the  finite  element  problem. 
Various  possible  implementations  of  the  method  used  to  recover  gradients 
are  possible  and  have  been  proposed  and  tested  in  the  past,  for  example 
the  superconvergent  patch  recovery  (SPR),  the  Clemant  interpolation,  and 
various  form  of  averaging  [19,  24]. 

Usually  these  procedures  amount  to  averaging  the  piecewice  constant  gra¬ 
dient  Vuh  over  suitable  patches  of  elements  around  a  given  node  (patch 
assembly  point),  in  order  to  obtain  a  continuous  piecewise  linear  recovered 
gradient.  This  recovered  value  is  sometimes  found  to  superconverge  to  the 
gradient  of  u  on  a  set  of  points,  which  coincide  with  the  patch  assembly 
point  when  using  linear  elements. 

The  concept  of  effectivity  index  is  frequently  used  to  measure  the  quality 
of  an  error  estimator.  The  effectivity  index  is  defined  as 


where  rj  is  the  global  error  indicator 
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over  the  triangulation  Th  of  the  domain  f l,  t]k  being  the  elemental  error 
indicator.  An  error  estimator  is  said  to  be  asymptotically  exact  if  0  ->  1 
when  the  mesh  size  h  approaches  zero.  There  is  a  relationship  between  the 
accuracy  of  the  recovered  gradients  and  the  effectivity  index,  and  in  fact  if 
the  gradient  recovery  is  superconvergent,  the  error  estimator  is  asymptoti¬ 
cally  exact  [29,  30].  The  ZZ  estimator  is  indeed  asymptotically  exact,  if  the 
solution  is  sufficiently  smooth  and  the  grid  possesses  certain  properties  of 
regularity  [24].  In  practice,  the  ZZ  estimator  has  been  applied  to  more  gen¬ 
eral  unstructured  grids  with  interesting  results  [30].  This  indeed  motivates 
the  use  of  the  ZZ  procedures  even  in  the  context  of  the  applications  that  are 
releveant  to  this  research  project. 

It  should  be  mentioned  that  in  reality,  for  realistic  adaptive  applications 
where  the  grid  size  will  clearly  not  tend  to  zero,  or  on  highly  anisotropic 
grids,  the  requirement  of  asymptotical  exactness  is  more  appropriately  re¬ 
placed  by  a  saturation  assumption  [27],  which  can  be  expressed  as 

II  Vzzuh  _  Vu  ||i2(n)<  P  ||  -  Vu  ||L2(n), 
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with  /3  e  [0, 1),  where  Vzzuh  is  the  ZZ  recovered  (improved)  gradient. 

In  this  work  we  use  the  ZZ  estimator  for  driving  the  adaption  process. 
At  first,  we  have  generalized  the  ZZ  method  to  the  compressible  Euler  and 
Navier-Stokes  equations  written  in  terms  of  entropy  variables.  This  yields 
a  numerical  procedure  that  is  directly  usable  per-se  and  that  has  been  im¬ 
plemented  in  our  parallel  adaptive  software  based  on  the  Galerkin  Least- 
Squares  finite  element  formulation  [9,  8].  This  formulation  of  the  estimator 
is  briefly  described  in  Section  3.2.  Next,  we  try  to  go  beyond  the  simple  use 
of  the  ZZ  estimator,  and  we  present  some  preliminary  work  on  the  coupling 
of  the  ZZ  recovery  with  an  estimate  in  the  energy  norm  of  the  interpolation 
error  which  takes  into  account  explicitely  the  anisotropic  information  of  the 
solution  and  of  the  mesh  elements.  This  part  of  the  work  is  described  is 
Section  3.3. 


3.2.  An  Implementation  of  the  ZZ  Estimator  for  the  Compressible 
Euler  and  Navier-Stokes  Equations.  The  Navier-Stokes  equations  can 


be  written  as 

(2)  A0 V,t  +  AiV,*  =  (K#V j)  .  +  F. 


Denoting  by  U  the  conservative  variables,  V  =  dH/dU  are  the  entropy 
variables,  defined  starting  form  a  generalized  entropy  function  H  =  H( U). 
Furthermore,  Ao  =  9U / 9W  denotes  the  associated  Reimannian  metric  ten- 
sor,  which  can  be  used  for  defining  consistent  energy  norms. 

To  discretize  the  Navier-Stokes  equations  in  space  and  time,  we  use  the 
Time-Discontinuous  Galerkin  finite  element  method,  stabilized  using  the 
least-squares  formulation  of  ref.  [26].  The  spatial  discretization  uses  linear 
shape  functions. 

According  to  the  ZZ  method,  the  energy  norm  of  the  error  is  then  com¬ 
puted  as 

(3)  ||  eh  ||=  (j  {vZZVh  ~  W'1)  •  diag[A0]  ( VzzVh  -  VVft)  dnj  , 


where 


diag[A0]  = 


Ao 


k-o  J 


and 


V(-) 


I 

dx\ 

aio 

.  dxd  - 


where  d  is  the  number  of  space  dimensions.  VzzVh  is  the  recovered  gradient 
of  the  entropy  variables,  which  is  here  computed  based  on  averaging  using 
element  volumes.  This  technique  is  here  favoured  to  the  alternative  offered 
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by  least-squares  fitting.  In  fact,  it  is  not  uncommon  in  realistic  applications 
to  incur  in  local  mesh  configurations  that  will  make  the  least-squares  re¬ 
construction  singular  or  very  ill-conditioned.  For  example,  this  will  happen 
every  time  one  has  the  centroids  of  the  patch  of  regions  that  lie  on  the  same 
plane. 

The  implementation  of  the  formulation  is  straightforward  and  was  easily 
parallelized  using  the  message  passing  protocol  MPI,  as  for  all  the  rest  of 
the  code. 

This  form  of  the  estimator  can  be  used  for  driving  isotropic  adaption  pro¬ 
cesses.  Elements  where  the  local  percentage  error  is  above  a  user-determined 
threshold  are  refined.  In  this  work  we  use  edge-based  mesh  modification 
primitives,  and  therefore  the  elemental  errors  are  reinterpolated  onto  the 
mesh  edges  before  adaption.  In  turn,  edges  marked  for  refinement  are  split 
using  all  possible  tetrahedral  edge-splitting  templates.  Alternatively,  edges 
can  be  collapsed  if  local  coarsening  is  requested. 

3.2.1.  Numerical  Examples.  We  report  here  two  simple  numerical  experi¬ 
ments  that  were  used,  among  others,  to  assess  the  correct  implementation 
of  the  methodology. 

The  first  example  is  a  two-dimensional  steady  problem  consisting  of  a 
Mach  two  flow  over  a  wedge  at  an  angle  of  10  deg,  resulting  in  the  creation 
of  an  oblique  shock  at  the  wedge  leading  edge.  Since  the  code  only  supports 
three-dimensional  elements,  the  problem  was  solved  adding  a  third  dimen¬ 
sion  equal  to  one  tenth  of  the  domain  size.  Figure  25  gives  the  true  (at  left) 
and  computed  (at  right)  error  on  the  initial  coarse  grid.  Good  matching 
between  the  two  quantities  can  be  observed  at  all  locations  in  the  regions 
of  the  shock.  Next,  an  automated  adaptive  analysis  was  started,  using  the 
computed  ZZ  error  to  drive  the  mesh  refinement  as  previously  explained. 
Figure  26  shows  the  obtained  grid  after  two  refinement  steps.  A  nice  and 
well  distributed  clustering  of  elements  can  be  observed  in  the  region  of  the 
shock. 

Next,  we  consider  the  classical  Onera  M6  wing  in  transonic  fligh.  The 
wing  is  characterized  by  an  aspect  ratio  of  3.8,  a  leading  edge  sweep  angle  of 
30  deg,  and  a  taper  ratio  of  0.56.  The  airfoil  section  is  an  Onera  D  symmetric 
section  with  10%  maximum  thickness-to-cord  ratio.  We  consider  a  steady 
flow  problem  characterized  by  an  angle  of  attack  a  =  3.06  deg  and  a  value  of 
M  =  0.8395  for  the  freestream  Mach  number.  In  such  conditions,  the  flow 
pattern  around  the  wing  is  characterized  by  a  complicated  double-lambda 
shock  on  the  upper  surface  of  the  wing  with  two  triple  points. 

Figure  27  shows  the  error  computed  by  the  implemented  error  estimator 
on  the  initial  grid.  The  error  correctly  identifies  the  region  of  the  shock,  and 
also  targets  the  leading  and  trailing  edges  for  refinement.  Finally,  figure 
28  shows  the  adapted  mesh  obtained  after  two  analysis-estimation-adaption 
cycles. 
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3.3.  An  Anisotropic  Recovery-Based  A-Posteriori  Error  Estima¬ 
tor.  In  this  section  we  report  an  alternative  use  of  the  ZZ  recovery  tech¬ 
nique,  that  aims  at  developing  a  true  anisotropic  estimator.  It  should  in 
fact  be  noticed  that  the  ZZ  method,  in  its  original  formulation,  does  not 
directly  allow  for  the  construction  of  a  local  metric  [16]  that  can  then  be 
used  for  driving  an  anisotropic  adaption  procedure.  The  basic  idea  of  the 
formulation  is  briefly  described  in  the  following,  and  more  details  can  be 
found  in  ref.  [22]  and  the  literature  cited  therein. 

At  first,  we  define  a  linear  mapping  that  associates  the  reference  element 
K  to  the  actual  element  K : 

(4)  x  =  Mkx  +  rh. 

The  non-singular  Jacobian  Mk  can  be  factored  by  a  polar  decomposition  as 

(5)  Mk  —  BkRk, 

where  Bk  is  symmetric  and  positive  definite,  while  Rk  is  orthogonal  and 
amounts  to  a  rotation.  The  eigenvalues  of  Bk  axe  denoted  \,k,  *  =  1?^ 
and  their  associated  eigenvectors  are  denoted  a^K,  i  =  l,d. 

It  is  possible  to  prove  under  appropriate  conditions  that,  given  the  bilinear 
form  B(uh,vh)  corresponding  to  the  weak  form  of,  e.g.,  a  diffusion  problem, 
the  bilinear  form  B(eh,vh),  where  eh  =  u  -  uh  is  again  the  discretization 
error,  can  be  bounded  as  follows 

(6)  \\B(e\v)\\  <  OtKPK{uh)wK{v). 

KeTh 

The  first  quantity  appearing  in  the  previous  expressions  can  be  synthetically 
given  as 

oik  =  ock{\k), 

while  the  second  term  is 

pK{uh)  =  PK(rK(uh),  RK(uh),  min(Ai>K)). 


In  the  previous  expression,  denotes  the  elemental  internal  residuals, 

while  RK(uh)  represents  the  elemental  boundary  residuals.  Finally,  the  last 
term  writes 

wK(v)  =  wK{GK{v),\i,K,ai,K), 
where  Gk(v )  is  a  symmetric  non-negative  matrix  defined  as 


Vt>  <g>  Vv  d T, 


where  Vh  is  the  patch  of  all  elements  that  share  a  vertex  with  element  K. 

These  results  can  be  used  for  computing  a-posteriori  error  estimates  based 
on  the  solution  of  the  dual  problem.  This  can  be  very  expensive  for  non  self- 
adjoint  problems.  A  more  interesting  idea  is  to  set  v  =  eh  in  (6),  to  get  an 
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implicit  estimate  for  the  square  of  the  energy  norm  of  the  discretization  error 
as 

(7)  ||  eh  ||<  ^2  aKpK(uh)wK(eh). 

Ken 

Since  the  right  hand  side  of  (7)  still  depends  on  the  error,  we  need  to  replace 
wK(eh)  with  a  computable  quantity.  An  elegant  solution  is  to  use  the  ZZ 
recovery  technique  to  replace  Gn(fi^1')  with  a  new  quantity  Gk  defined 

as 

GK{ezz’h )  =  Y,  /  &ZZvh  ~  0  (V ZZyh  ~  Vv^ dT’ 

Tevh Jt 

Finally  we  get  a  computable  a-posteriori  recovery-based  error  estimator  de¬ 
fined  on  each  mesh  element  as 

(8)  tj%  =  ^2  <XKPK{uh)wK{ezz’h). 

Ken 

This  estimator  provides  anisotropic  information  in  a  consistent  manner, 
since  it  is  based  on  the  ZZ  recovery  that  is  known  to  perform  very  well 
in  a  number  of  circumstances,  and  is  also  cheap  to  compute. 

Ref.  [22]  gives  details  on  how  it  is  possible  to  reconstruct  a  local  metric 
based  on  the  estimator,  in  order  to  drive  an  anisotropic  mesh  modification 
procedure.  In  particular,  a  piecewise  constant  metric  matrix  can  be  com¬ 
puted  on  each  element  of  the  triangulation  based  on  the  knowledge  of  its 
eigenvalues  and  eigenvectors.  These  quantities  can  be  computed  based  on 
the  estimator  (8),  provided  certain  requirements  on  the  error  distribution 
and  on  the  grid  are  specified.  In  particular,  it  makes  sense  to  request  that 
the  error  be  equidistributed  throughout  the  grid  and  that  the  number  of 
elements  is  minimized.  This  leads  to  a  minimization  problem  that  gives  the 
required  definition  of  the  metric  matrices. 

Preliminary  results  concerning  the  proposed  technique  are  presented  in 
ref.  [22]  on  academic  test  cases. 
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4.  Conclusions 

In  this  document,  we  have  summarized  the  research  conducted  to  date  on 
anisotropic  adaption  using  unstructured  finite  elements,  under  the  support 
of  the  European  Research  Office  of  the  U.S.  Army.  This  work  was  focused  in 
two  main  areas:  the  development  of  anisotropic  mesh  adaption  procedures 
and  the  problem  of  a-posteriori  error  estimation.  These  two  points  are  also 
reflected  in  the  two  tasks  (1.1  and  1.2)  of  the  present  contract,  which  read. 

•  Task  1.1:  Analysis  and  development  of  anisotropic  error  estimation 
procedures. 

•  Task  1.2:  Analysis  and  development  of  anisotropic  mesh  refinement 
strategies. 

The  delivarable  for  this  part  of  the  work  stated: 

Delivarable  (12  month  milestone):  technical  document  re¬ 
porting  on  the  results  of  the  research  in  the  field  of  anisotropic 
error  estimation  and  mesh  adaption,  suitable  for  rotorcraft 
problems.  The  document  will  also  detail  the  results  gath¬ 
ered  during  preliminary  testing  and  implementation  of  the 
methods,  with  the  goal  of  determining  their  suitability  for 
the  application  to  full  scale  hovering  rotorcraft  problems. 

Regarding  the  first  area  of  investigation,  we  have  here  presented  a  pro¬ 
cedure  for  the  insertion  of  anisotropic  tetrahedral  grid  layers  in  the  interior 
and  in  the  proximity  of  the  boundary  of  curved  three-dimensional  domains. 
The  present  approach  has  some  strenghts  and  some  weaknesses,  as  probably 
most  methods  do. 

Among  the  strong  points,  we  can  mention  here  that  it  allows  to  generate 
locally  adapted  anisotropic  meshes  with  excellent  control  of  the  grid  quality. 
In  particular,  since  the  anisotropic  layers  are  “almost”  structured  in  nature, 
the  stretched  elements  are  denoted  by  dihedral  angles  which  are  either  very 
small  or  very  close  to  ninety  degrees.  Theoretical  results  show  that  large 
dihedral  angles  affect  the  discretization  error  and  the  conditioning  of  the 
discrete  problem,  and  should  therefore  be  avoided.  It  should  be  noted  that 
other  anisotropic  refinement  procedures,  such  as  some  form  of  point  insertion 
based  on  local  anisotropic  metrics  [5],  can  not  guarantee  the  control  of  large 
dihedral  angles  and  in  fact  in  general  will  produce  very  poor  grids  with 
respect  to  this  quality  measure. 

The  present  method  allows  to  use  any  existing  mesh  generator  for  produc¬ 
ing  the  isotropic  grid,  and  is  not  limited  to  advancing  front  methods  as  in 
the  advancing  layers  procedure.  Furthermore,  the  method  does  not  require 
a  complex  post-processing  to  check  and  fix  cross-overs  and  self-intersections. 
The  software  procedures  are  integrated  with  a  solid  modeler,  that  provides 
the  required  geometric  and  topological  information. 

On  the  other  hand,  it  is  also  clear  that,  by  its  very  design,  the  present 
approach  is  only  suited  to  the  refinement  of  one-dimensional  anisotropic 
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features,  i.e.  for  solution  features  that  are  characterized  by  a  single  local  di¬ 
rection  of  high  gradients,  as  it  is  the  case  of  boundary  layers  and  shocks.  The 
generalization  of  the  ideas  here  reported  to  also  support  the  refinement  of 
two-dimensional  features,  such  as  vortices,  leads  to  very  complicated  mesh¬ 
ing  processes  which  are  probably  difficult  to  justify  for  practical  applications. 
In  this  sense,  it  seems  that  a  better  approach  would  be  to  complement  the 
present  procedures  with  a  metric-based  insertion  method  to  be  used  when 
and  where  the  layer  insertion  method  here  proposed  can  not  be  profitably 
used.  We  are  indeed  already  working  in  this  direction,  in  order  to  expand 
the  applicability  of  the  present  software  to  handle  these  more  general  cases. 

The  method  here  proposed  seems  well  suited  for  the  analysis  of  hovering 
rotors,  with  the  limitations  noted  above.  The  method  is  ideally  suited  to 
the  growth  of  boundary  layers  for  Navier-Stokes  simulations,  both  around 
the  blades  and  the  fuselage.  For  the  refinement  of  certain  internal  solution 
features,  as  for  example  in  the  analysis  of  the  acoustic  pressure  field  for 
hovering  rotors,  as  shown  in  figure  29,  the  region  of  localization  is  also  well 
suited  to  the  application  of  the  present  procedure.  However,  we  remark 
once  again  that  a  greater  flexibility  would  be  obtained  by  adding  to  this 
existing  capability  the  possibility  of  refining  by  metric-based  point  insertion, 
for  example  in  the  wake  and  tip-vortex  regions.  Future  work  will  in  fact 
concentrate  the  research  in  this  specific  direction. 

Regarding  the  area  of  a-posteriori  error  estimation,  we  have  reported  on 
a  generalization  of  the  ZZ  patch-recovery  technique  to  the  case  of  the  Euler 
and  Navier-Stokes  equations.  This  method  performs  extremely  well  on  reg¬ 
ular  meshes,  but  can  also  be  used  with  success  on  the  irregular  meshes  that 
characterize  realistic  unstructured  applications.  Furthermore,  the  method  in 
its  present  implementation  is  very  robust,  cheap  to  compute  and  paralleliz- 
able  with  great  efficiency  being  local  in  nature.  In  this  sense,  this  method 
seems  to  offer  a  mature  methodology  for  the  estimation  of  the  error  asso¬ 
ciated  with  finite  element  solutions  of  complex,  large  scale  problems.  This 
allows  to  avoid  the  use  of  gradient  norms  for  driving  the  adaption  process, 
that,  providing  no  information  on  the  finite  element  error  but  simply  track¬ 
ing  regions  of  high  solution  variation,  can  not  be  fully  automated.  We  have 
tested  with  success  the  procedure  on  a  number  of  examples,  and  we  were 
able  to  include  it  in  the  parallel  adaptive  code  with  no  particular  difficulty. 
Furthermore,  is  has  been  experimentally  observed  that  the  method  behaves 
well  even  in  the  presence  of  discontinuities,  which  is  the  typical  case  of  Euler 
shocks.  While  the  mathematical  theory  behind  this  recovery-based  estima¬ 
tor  will  clearly  not  hold  in  this  case,  the  procedure  is  still  able  to  identify 
the  regions  that  need  to  be  better  resolved. 

Although  the  ZZ  method  has  a  number  of  useful  characteristics,  it  was 
here  noted  that  it  does  not  provide  per-se  information  for  constructing  metric 
fields  defined  on  the  computational  domain  to  be  used  for  driving  anisotropic 
adaption  processes.  To  overcome  this  limitation,  we  have  proposed  to  resort 
to  interpolation  estimates  that  include  the  information  on  the  anisotropy  of 
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the  solution  as  well  as  the  anisotropy  of  the  computational  grid.  With  the  use 
of  such  estimates,  one  obtains  a  definition  of  a  local  metric  that  can  be  used 
in  an  iterative  fashion  to  obtain  an  anisotropically  adapted  locally  optimal 
mesh.  Unfortunately,  such  estimates  are  prohibitively  expensive  to  compute 
for  realistic  non-linear  non-self  adjoint  large  scale  applications,  like  the  ones 
that  are  relevant  to  this  research  project.  To  overcome  this  difficulty,  we 
have  proposed  to  exploit  the  sinergy  between  the  anisotropic  estimates  and 
the  ZZ  recovery  technique.  This  way,  one  can  compute  an  approximation 
to  the  local  error  appearing  in  the  anisotropic  estimates  using  the  robust 
and  cheap  ZZ  method,  making  the  anisotropic  information  readily  available 
for  use  in  a  fully  automated  adaption  process.  So  far  we  have  analyzed  this 
approach  from  a  mathematical  point  of  view  and  we  have  numerically  tested 
it  on  academic  examples.  The  results  so  obtained  are  very  promising  and 
encourage  to  continue  the  work  in  this  direction  towards  the  generalization 
of  these  ideas  to  the  compressible  Euler  and  Navier-Stokes  equations.  At 
present  we  do  no  have  results  to  corroborate  these  preliminary  findings,  but 
we  intend  to  pursue  the  research  in  this  direction  in  order  to  clarify  the 
suitability  of  the  aproach  for  applications  of  engineering  interest. 

Comparing  the  proposed  tasks  and  milestone  with  the  achieved  results, 
we  can  say  that  we  have  analyzed  and  tested  the  methods  here  proposed 
on  meaningful  examples,  and  we  have  already  implemented  all  the  most 
mature  procedures  in  the  parallel-adaptive  rotor  code.  This  part  of  the  work 
would  have  been  tasks  2.1  and  2.2  of  a  possible  second  year  continuation 
of  this  contract,  according  to  the  original  proposal.  Parts  which  are  less 
mature,  and  in  particular  the  anisotropic  estimators,  are  beeing  analyzed 
using  simpler  research  codes. 

The  results  of  this  research  activity  have  been  reported  in  the  following 
scientific  publications: 

•  Bottasso,  C.L.,  Detomi,  D.  (2002)  A  procedure  for  tetrahedral  bound¬ 
ary  layer  mesh  generation.  Engineering  with  Computers,  accepted, 
to  appear. 

•  Micheletti,  S.,  Perotto,  S.  (2001)  An  anisotropic  recovery-based  a- 
posteriori  error  estimator.  Proceedings  of  ENUMATH  2001,  Euro¬ 
pean  Conference  on  Numerical  Mathematics  and  Advanced  Appli¬ 
cations. 

•  Bottasso,  C.L.,  Detomi,  D.  (2002)  Procedures  for  the  refinement 
and  generation  of  anisotropic  tetrahedral  grids  in  three-dimensional 
domains.  SIMAI  2002,  Chia,  CA,  Italy,  May  27-31,  2002. 

•  Bottasso,  C.L.,  Detomi,  D.  (2002)  Tetrahedral  boundary  layer  mesh 
generation  for  viscous  flows  in  complex  geometries.  XVI  Congresso 
Nazionale  AIDAA,  Palermo,  Italy,  September  24-28,  2001. 

•  Bottasso,  C.L.  Detomi,  D.  (2002)  Generation  and  refinement  of 
anisotropic  tetrahedral  and  hybrid  grids  in  three  dimensions.  ANIS- 
GRID  2002,  Politecnico  di  Milano,  Italy,  June  21,  2002. 
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•  Formaggia,  L.,  Micheletti,  S.  (2002)  An  anisotropic  framework  for 
error  estimates  in  FEM.  ANISGRID  2002,  Politecnico  di  Milano, 
Italy,  June  21,  2002. 

•  Bottasso,  C.L.,  Detomi,  D.  (2002)  Insertion  of  anisotropic  layers  in 
three-dimensional  tetrahedral  grids.  In  preparation  for  submission. 

Finally,  we  remark  that  in  June  2002  we  have  organized  a  one-day  work¬ 
shop  at  the  Politecnico  di  Milano  on  the  topic  of  anisotropic  grids: 
ANISGRID  2002.  Anisotropic  Grids:  Generation,  Adap¬ 
tion  and  Error  Estimation.  Organizers:  C.L.  Bottasso,  S. 
Micheletti,  S.  Perotto,  R.  Sacco.  Politecnico  di  Milano,  Mi¬ 
lano,  Italy,  June  21,  2002. 

Speakers  and  participants  were  from  the  European  academic  and  industrial 
worlds.  This  initiative,  although  not  directly  financed  by  this  research  con¬ 
tract,  was  made  possible  by  the  research  activity  described  herein  and  was 
directly  inspired  by  it. 
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Figure  2.  Anisotropic  adaption  by  mesh  cutting,  opening 
and  filling.  Top:  original  mesh  and  trace  of  the  surface.  Mid¬ 
dle:  mesh  after  the  cut  and  the  local  optimization.  Bottom: 
final  adapted  grid. 
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Figure  3.  Automatic  curving  of  growth  curves  through 
mesh  relaxation  in  the  boundary  layer  region. 


Figure  4.  Growth  of  a  boundary  layer  mesh  in  a  sharp  cor¬ 
ner  region. 


Figure  7.  Additional  linear  springs  connecting  each  vertex 
with  the  opposite  face  in  a  tetrahedron. 
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Figure  13.  Internal  layer  ending  within  the  domain,  i.e. 
without  intersecting  a  model  boundary. 
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FIGURE  14.  Boundary  layer  spillover,  a)  Required  growth 
directions  without  spillover;  b)  patch  of  mesh  faces  that  will 
be  detached  from  the  model;  c)  detachment  and  creation  of 
boundary  layer  with  spillover. 
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Figure  15.  Classification  against  model  entities  of  the 
newly  created  mesh  entities. 


FIGURE  16.  Boundary  layer  mesh  grown  around  an  engine 
nacelle  model.  Parts  of  the  anisotropic  grid  were  deleted 
for  exposing  the  through-the-thickness  discretization  and  the 
original  model  face  triangulation. 


<85  85-95  95-105  105-115  115-125  125-135  135-145  145-155  155-165  >165 


FIGURE  17.  Histogram  of  the  maximum  dihedral  angles  for 
each  element  of  the  initial  (thin  lines)  and  final  (thick  lines) 
meshes  of  the  nacelle  model. 


Figure  20.  Boundary  layer  mesh  grown  around  a  submarine 
model.  Parts  of  the  anisotropic  grid  were  deleted  for  expos¬ 
ing  the  through-the-thickness  discretization  and  the  original 
model  face  triangulation. 


Figure  21.  Details  of  a  boundary  layer  mesh  grown  around 
a  submarine  model. 
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Figure  23.  Insertion  of  a  thin  internal  layer  in  front  of  a  sphere. 
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FIGURE  24.  Insertion  of  a  thin  internal  layer  in  front  of  a 
sphere.  Elements  in  the  isotropic  grid  were  stripped  away  in 
order  to  show  the  anisotropic  mesh. 
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Figure  25.  True  (left)  and  ZZ  estimated  (right)  errors  for 
the  oblique  shock  problem. 
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Figure  26.  Adapted  grid  after  two  refinements  for  the 
oblique  shock  problem. 
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Figure  27.  ZZ  estimated  error  for  the  Onera  M6  wing  problem. 
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Figure  28.  Adapted  grid  after  two  refinements  for  the  On- 
era  M6  wing  problem. 
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Figure  29.  View  of  an  isotropically  refined  mesh  for  a  high 
speed  impusive  noise  computation. 


